// ContinuousTabuSearch implementation
/////////////////////////////////////////////////////////////////
#include <iomanip>
#include <vector>

#include "continuous_tabusearch.h"
#include "hypersphere.h"
#include "tabulist.h"
#include "wrapperFunctions.h"
#include "molecule.h"

#include <algorithm>

#define PI   3.1415926535897932385



//////////////////////////////////////////////////////////////////
// Continuous Tabu Search : constructors/destructors
//////////////////////////////////////////////////////////////////

ContinuousTabuSearch::ContinuousTabuSearch(int numDimension, int numNeighSol, double maxRadius, double tabuRadius,
                                           double promRadius, int maxTabuListSize, int maxPromlistSize,
                                           Function *Func) : 
  m_Tabulist(maxTabuListSize), m_Promlist(maxPromlistSize), m_neighSolution(numNeighSol, vector<double> (numDimension))
{
  // initialise parameters
  m_NumDimension = numDimension;
  m_NumNeighSol = numNeighSol;
  m_maxTabulistSize = maxTabuListSize;
  m_maxPromlistSize = maxPromlistSize;

  // Used to convert from spherical to cartesian coords
  newA(m_Phi, m_NumDimension - 1);
  newA(m_Cos, m_NumDimension);
  newA(m_Ran, m_NumDimension);
  m_Sin = new double[m_NumDimension];

  // the function to be used
  StoredFunction = Func;  
}

void ContinuousTabuSearch::Initialise(double maxRadius, double tabuRadius, double promRadius)
{
  double *solution;
  newA(solution, m_NumDimension);

  m_maxRadius = maxRadius;
  m_promRadius = promRadius;
  m_tabuRadius = tabuRadius;

  nFuncNoCall = 0;
  nPromNoCall = 0;
  nReduction = 0;

  m_Promlist.Reset(m_maxPromlistSize);
  m_Tabulist.clear(m_maxTabulistSize);

  for(int i = 1; i<=m_maxPromlistSize; i++)
  {
    StoredFunction->createRandomSol(solution);
    funcval = StoredFunction->Value(solution);
    //funcval = StoredFunction->NormValue(solution);
    
    if( m_Promlist.isInList(solution) )
    {
      i--;
    }
    else
    {
      m_Promlist.Initialise(m_NumDimension, m_maxPromlistSize, solution, m_promRadius, funcval);
    }
  }

  deleteA(solution);

  createCube(m_maxRadius, m_tabuRadius);
  //createRadii(m_maxRadius, m_tabuRadius);
}

//////////////////////////////////////////////////////////////////
// Hybrid Monte Carlo : functions
//////////////////////////////////////////////////////////////////

void ContinuousTabuSearch::doRun(double *solution, int maxIterations, string& directory)
{
  for (int i = 1; i <= maxIterations; i++)
  {
    thresholdValue = m_Promlist.Threshold(m_maxPromlistSize);

    // Create a set of neighbouring solutions
    neighbour(solution);
    origValue = StoredFunction->Value(solution);
    //origValue = StoredFunction->NormValue(solution);    // gets the function value of

    //origValue = StoredFunction->Result.getBestValue();
    StoredFunction->nCall--;

    // Sets the best neighbouring solution to be the first one
    bestForThisIteration = StoredFunction->Value(m_neighSolution[0]);
    //bestForThisIteration = StoredFunction->NormValue(m_neighSolution[0]);
    worstForThisIteration = bestForThisIteration;
    int bestNeighSol = 0;   // Keeps track of the neighbouring
                            // solution which gives the best value
    int worstNeighSol = 0;

    double aValue;

    // loop over the remaining solutions to see which is the best
    // and store the value and which neighbouring solution gives
    // this solution
    for (int j = 1; j < m_NumNeighSol; j++)
    {
      aValue = StoredFunction->Value(m_neighSolution[j]);
      //aValue = StoredFunction->NormValue(m_neighSolution[j]);
      if (bestForThisIteration > aValue)
      {
        bestForThisIteration = aValue;
        bestNeighSol = j;
      }
      else if (worstForThisIteration < aValue)
      {
        worstForThisIteration = aValue;
        worstNeighSol = j;
      }
    }

    m_Tabulist.add(m_NumDimension, solution, m_tabuRadius);

    // if f(s1)=bestForThisIteration < f(s2) < ... < f(sn)
    // then if f(s) < f(s1) < f(s2) < ... < f(sn)
    // then invoke the simplex/CGM
    if (bestForThisIteration > origValue)
    {
      nFuncNoCall++;

      if ( bestForThisIteration > thresholdValue )
      {
        nPromNoCall = 0;
        nReduction++;
        //cgmBestForThisIteration = bestForThisIteration;
        simpBestForThisIteration = bestForThisIteration;
//        wrapperCGM(solution,m_NumDimension,1e-3,&cgmBestForThisIteration,StoredFunction, directory);
        wrapperSimplex(solution, m_NumDimension, 1e-3, 1e-2, simpBestForThisIteration, StoredFunction);

        // add to promising list and clear tabu list
        m_Promlist.replace(m_NumDimension, solution, m_promRadius, simpBestForThisIteration);
        m_Tabulist.clear(m_maxTabulistSize);
        m_maxRadius /= 2;
        m_tabuRadius /= 2;
        m_promRadius /= 2;
        createCube(m_maxRadius, m_tabuRadius);
        //createRadii(m_maxRadius, m_tabuRadius);
      }
      else
      {
        nPromNoCall++;
        nReduction = 0;
        
        // update the new solution with the best neighbouring solution
        for (int k = 0; k < m_NumDimension; k++)
        {
          solution[k+1] = m_neighSolution[bestNeighSol][k];
        }
      }
    }
    else
    {
      nPromNoCall++;
      nFuncNoCall = 0;
      nReduction = 0;
      
      // update the new solution with the best neighbouring solution
      for (int k = 0; k < m_NumDimension; k++)
      {
        solution[k+1] = m_neighSolution[bestNeighSol][k];
      }
    }

    /*if(nPromNoCall > 2*m_NumDimension)
    {
      break;
    }

    if(nFuncNoCall > 5*m_NumDimension)
    {
      break;
    }

    if(nReduction > 2*m_NumDimension)
    {
      break;
    }*/
  }
}

void ContinuousTabuSearch::doRunDiagnosis(double *solution, int maxIterations, char filename[], char filecoord[], string& directory)
{
  ofstream outputcoords(filecoord);
  outputcoords.precision(6);
  outputcoords << setw(4) << 0 << setw(13);

  for(int l = 1; l<=m_NumDimension; l++)
  {
    outputcoords << solution[l] << setw(14);
  }
  outputcoords << endl;

  // Open diagnostic file and write the names of the columns
  ofstream funcfile(filename);
  funcfile.precision(6);
  funcfile << "Iteration | Start Func Value |  Best Neighbour Value  |  Func Value  " 
    << "|  Function Evaluations  |  Distance from known Global"
    << endl << setw(11) << "|" << setw(14) << "(bef min)" << setw(5) << "|"
    << setw(25) << "|"
    << setw(11) << "(aft min)" << setw(4) << "|" << setw(25) << "|" << setw(18) << "Optimum" << endl
    << "------------------------------------------------------------------------------------------"
    << "----------------------------------" << endl
    << setw(11) << "|" << setw(19) << "|" << setw(25) << "|" << setw(15) << "|" << setw(25) << "|"
    << endl;
  
  for (int i = 1; i <= maxIterations; i++)
  {
    outputcoords << setw(4) << i << setw(13);

    // Create a set of neighbouring solutions
    neighbour(solution);
    origValue = StoredFunction->Value(solution);
    //origValue = StoredFunction->NormValue(solution);    // gets the function value of

    thresholdValue = m_Promlist.Threshold(m_maxPromlistSize);//, origValue);

    //origValue = StoredFunction->Result.getBestValue();
    StoredFunction->nCall--;
    
    funcfile << setw(6) << i 
      << setw(17) << origValue;

    // Sets the best neighbouring solution to be the first one
    bestForThisIteration = StoredFunction->Value(m_neighSolution[0]);
    //bestForThisIteration = StoredFunction->NormValue(m_neighSolution[0]);
    worstForThisIteration = bestForThisIteration;
    int bestNeighSol = 0;   // Keeps track of the neighbouring
                            // solution which gives the best value
    int worstNeighSol = 0;

    double aValue;

    // loop over the remaining solutions to see which is the best
    // and store the value and which neighbouring solution gives
    // this solution
    for (int j = 1; j < m_NumNeighSol; j++)
    {
      aValue = StoredFunction->Value(m_neighSolution[j]);
      //aValue = StoredFunction->NormValue(m_neighSolution[j]);
      if (bestForThisIteration > aValue)
      {
        bestForThisIteration = aValue;
        bestNeighSol = j;
      }
      else if (worstForThisIteration < aValue)
      {
        worstForThisIteration = aValue;
        worstNeighSol = j;
      }
    }

    funcfile << setw(22) << bestForThisIteration;
    m_Tabulist.add(m_NumDimension, solution, m_tabuRadius);

    // if f(s1)=bestForThisIteration < f(s2) < ... < f(sn)
    // then if f(s) < f(s1) < f(s2) < ... < f(sn)
    // then invoke the simplex/CGM
    if (bestForThisIteration > origValue)
    {
      nFuncNoCall++;
      double sumneigh = 0.0;

      for(int t = 0; t < m_NumNeighSol; t++)
      {
        sumneigh += StoredFunction->Value(m_neighSolution[t]);

        StoredFunction->nCall--;
      }
      //sumneigh /= m_NumNeighSol;
      //sumneigh -= origValue;

      if ( sumneigh > fabs(thresholdValue-origValue) )
      {
        nPromNoCall = 0;
        nReduction++;
        //doshSimpBestForThisIteration = bestForThisIteration;
        cgmBestForThisIteration = bestForThisIteration;
        //simpBestForThisIteration = bestForThisIteration;
        wrapperCGM(solution,m_NumDimension,1e-3,&cgmBestForThisIteration,StoredFunction, directory);
        //wrapperSimplex(solution, m_NumDimension, 1e-3, 1e-2, simpBestForThisIteration, StoredFunction);
        //wrapperDoshSimplex(solution, 0.1, m_NumDimension, doshSimpBestForThisIteration, StoredFunction); 

        funcfile << setw(20) << cgmBestForThisIteration << setw(17)
          << StoredFunction->nCall << setw(30) << StoredFunction->distanceFromGlobalMinimum(solution) << endl;

        for (int k = 0; k < m_NumDimension; k++)
        {
          outputcoords << solution[k+1] << setw(14);
        }
        outputcoords << endl;

        // add to promising list and clear tabu list
        m_Promlist.replace(m_NumDimension, solution, m_promRadius, simpBestForThisIteration);
        m_Tabulist.clear(m_maxTabulistSize);
        //m_maxRadius  /= 2;
        //m_tabuRadius /= 2;
        //m_promRadius /= 2; 
        //createCube(m_maxRadius, m_tabuRadius);
        //createRadii(m_maxRadius, m_tabuRadius);
      }
      else
      {
        nPromNoCall++;
        nReduction = 0;

        funcfile << setw(37) << StoredFunction->nCall << setw(30)
          << StoredFunction->distanceFromGlobalMinimum(solution) << endl;
        
        // update the new solution with the best neighbouring solution
        for (int k = 0; k < m_NumDimension; k++)
        {
          solution[k+1] = m_neighSolution[bestNeighSol][k];
          outputcoords << solution[k+1] << setw(14);
        }
        outputcoords << endl;
      }
    }
    else
    {
      nPromNoCall++;
      nFuncNoCall = 0;
      nReduction = 0;

      funcfile << setw(37) << StoredFunction->nCall << setw(30)
        << StoredFunction->distanceFromGlobalMinimum(solution) << endl;
      
      // update the new solution with the best neighbouring solution
      for (int k = 0; k < m_NumDimension; k++)
      {
        solution[k+1] = m_neighSolution[bestNeighSol][k];
        outputcoords << solution[k+1] << setw(14);
      }
      outputcoords << endl;
    }

    /*if(nPromNoCall > 2*m_NumDimension)
    {
      break;
    }

    if(nFuncNoCall > 5*m_NumDimension)
    {
      break;
    }

    if(nReduction > 2*m_NumDimension)
    {
      break;
    }*/
  }

  funcfile.close();
  outputcoords.close();
}

// this is just a wrapper function to call tracker method best()
// to save the best function value and coordinates
double ContinuousTabuSearch::BestSoFar(double *solution)
{
  double x;

  StoredFunction->Result.Best(x, solution);

  return (x);
}

/*************************************************************************************/
// partition the hypersphere by creating smaller radii
void ContinuousTabuSearch::createRadii(double maxRadius, double minRadius)
{
  m_radii.clear();
  m_radii.resize(0);

  double partition = (maxRadius-minRadius)/m_NumNeighSol;

  // add the radii to m_radii with the largest radius first (maxRadius)
  // and the smallest radius last (minRadius)
  for (int i = 0; i <= m_NumNeighSol; i++)
  {
    //m_radii.push_back( (maxRadius*(m_NumNeighSol - i + 1)) / m_NumNeighSol );
    m_radii.push_back( (maxRadius - (i * partition)) );
  }
  //m_radii.push_back(minRadius);

  if ( m_radii.size() != m_NumNeighSol + 1 ) 
  {
    cout << "m_radii attribute in ContinuousTabuSearch does not have the right size\n";
    exit(1);
  }
}
/*************************************************************************************/

/*************************************************************************************/
// used in the hypercube method
// partition the hypercube by creating smaller cubes
void ContinuousTabuSearch::createCube(double maxEdge, double minEdge)
{
  m_radii.clear();
  m_radii.resize(0);

  double partition = (maxEdge-minEdge)/m_NumNeighSol;

  // add the radii to m_radii with the largest radius first (maxRadius)
  // and the smallest radius last (minRadius)
  for (int i = 0; i <= m_NumNeighSol; i++)
  {
    //m_radii.push_back( (maxEdge*(m_NumNeighSol - i + 1)) / m_NumNeighSol );
    m_radii.push_back( (maxEdge - (i * partition)) );
  }
  //m_radii.push_back(minEdge);

  if ( m_radii.size() != m_NumNeighSol + 1 ) 
  {
    cout << "m_radii attribute in ContinuousTabuSearch does not have the right size\n";
    exit(1);
  }
}
/*************************************************************************************/

// create neighbouring solutions by populating m_neighSolution
void ContinuousTabuSearch::neighbour(double *solution)
{
  // aim to create a neighbouring solution for each shell one at a time
  for (int i = 0; i < m_NumNeighSol; i++)
  {
    // pick a random point in a shell
    pickRandomPointInCube(solution, m_radii[i+1], m_radii[i], i);
    //pickRandomPointInShell(solution, m_radii[i+1], m_radii[i], i);

    // checks if this random point is tabu
    if( m_Tabulist.isInList(m_neighSolution[i]) || m_Promlist.isInList(m_neighSolution[i]) )
    {
      i--; 
    }
  }
}

/*************************************************************************************/
// picks a random point in a shell between innerRadius and outerRadius
void ContinuousTabuSearch::pickRandomPointInShell(double *solution, double innerRadius,
                                                  double outerRadius, int neighSol)
{
  // generates a random number between the innerRadius and outerRadius
  m_RanRadii = ran1()*(outerRadius - innerRadius) + innerRadius;

  // generates the random angles, where the first numDim-2 angles
  // are between 0 and Pi and the last numDim-1 angle is between
  // 0 and 2Pi
  for (int i = 1; i < (m_NumDimension-1); i++)
  {
    m_Phi[i] = ran1()*PI;
  }
  m_Phi[m_NumDimension-1] = ran1()*(2*PI);

  // convert the spherical point to a cartesian point
  convertCartesian(m_RanRadii, m_Phi, m_neighSolution[neighSol]);

  // adding the solution coordinates to each neighbouring coordinate
  // so that we obtain a neighbouring solution inside a spherical shell
  // with the solution at its centre.
  for (int j = 0; j < m_NumDimension; j++)
  {
    m_neighSolution[neighSol][j] += solution[j+1];
  }
}
/*************************************************************************************/

/*************************************************************************************/
// picks a random point in a shell between innerRadius and outerRadius
void ContinuousTabuSearch::pickRandomPointInCube(double *solution, double innerEdge,
                                                  double outerEdge, int neighSol)
{
  for (int i = 1; i <= m_NumDimension; i++)
  {
    m_Ran[i] = (2*ran1()) - 1;
  }

  double biggestRan = fabs(m_Ran[1]);
  for (int j = 2; j <= m_NumDimension; j++)
  {
    if(fabs(m_Ran[j]) > biggestRan)
    {
      biggestRan = fabs(m_Ran[j]);
    }
  }

  double m_Big = innerEdge / biggestRan;
  double m_Root = pow(m_Big, (double) m_NumDimension);
  m_Root += pow(outerEdge, (double) m_NumDimension);
  m_Root -= pow(innerEdge, (double) m_NumDimension);
  double root = 1.0/((double) m_NumDimension);
  double new_Root = pow(m_Root, root);
  for (int k = 0; k < m_NumDimension; k++)
  {
    m_neighSolution[neighSol][k] = new_Root * m_Ran[k+1];
  }

  // adding the solution coordinates to each neighbouring coordinate
  // so that we obtain a neighbouring solution inside a spherical shell
  // with the solution at its centre.
  for (int l = 0; l < m_NumDimension; l++)
  {
    m_neighSolution[neighSol][l] += solution[l+1];
  }
}
/*************************************************************************************/

/*************************************************************************************/
// convert a spherical point to a cartesian point
void ContinuousTabuSearch::convertCartesian(double randomRadii, double *angle, vector<double> &neighSolArray)
{
  // initialises all of the neighbouring points to zero
  for (int i = 0; i < m_NumDimension; i++)
  {
    neighSolArray[i] = 0.0;
  }

  // if the dimension is 1 then the neighbouring solution is just the randomRadii
  if(m_NumDimension == 1)
  {
    neighSolArray[0] = randomRadii;
  }

  // m_Sin is a standard array of size m_NumDimension
  // m_Cos is a numerical receipe array of size m_NumDimension
  m_Sin[0] = 1.0;
  for (int j = 1; j < m_NumDimension; j++)
  {
    m_Cos[j] = cos(angle[j]);
    m_Sin[j] = sin(angle[j]);
  }
  m_Cos[m_NumDimension] = 1.0;

  // converting the spherical coordinate to a cartesian coordinate
  neighSolArray[0] = randomRadii * m_Sin[0];
  for (int k = 1; k < m_NumDimension; k++)
  {
    neighSolArray[k] = neighSolArray[k-1] * m_Sin[k];
  }

  for (int l = 0; l < m_NumDimension; l++)
  {
    neighSolArray[l] *= m_Cos[l+1];
  }
}
/*************************************************************************************/

void ContinuousTabuSearch::testHypercube(double *solution, double innerEdge,
                                         double outerEdge, int neighSol)
{
  for (int i = 1; i <= m_NumDimension; i++)
  {
    m_Ran[i] = (2*ran1()) - 1;
  }

  for (int r =  1; r <= m_NumDimension; r++)
  {
    cout << "m_Ran[" << r << "] = " << m_Ran[r] << endl;
  }

  double biggestRan = m_Ran[1];
  for (int j = 2; j <= m_NumDimension; j++)
  {
    if(m_Ran[j] > biggestRan)
    {
      biggestRan = m_Ran[j];
    }
  }

  double m_Big = innerEdge / biggestRan;
  double m_Root = pow(m_Big, (double) m_NumDimension);
  m_Root += pow(outerEdge, (double) m_NumDimension);
  m_Root -= pow(innerEdge, (double) m_NumDimension);

  double root = 1.0/(double) m_NumDimension;
  double new_Root = pow(m_Root, root);
  for (int k = 0; k < m_NumDimension; k++)
  {
    m_neighSolution[neighSol][k] = new_Root * m_Ran[k+1];
  }

  // adding the solution coordinates to each neighbouring coordinate
  // so that we obtain a neighbouring solution inside a spherical shell
  // with the solution at its centre.
  for (int l = 0; l < m_NumDimension; l++)
  {
    m_neighSolution[neighSol][l] += solution[l+1];
  }

  for (int m = 0; m < m_NumDimension; m++)
  {
    cout << "m_neighSolution[" << neighSol << "][" << m << "] = " << m_neighSolution[neighSol][m] << endl;
  }

  for (int n = 1; n <= m_NumDimension; n++)
  {
    cout << "solution[" << n << "] = " << solution[n] << endl;
  }
}


//////////////////////////////////////////////////////////////////
// ContinuousTabuSearch : friends
//////////////////////////////////////////////////////////////////

ostream &operator<<(ostream &stream, ContinuousTabuSearch &ContTabu)
{
  for (int i = 1; i <= ContTabu.m_NumNeighSol; i++)
  {
    for (int j = 1; j <= ContTabu.m_NumDimension; j++)
    {
      stream << "m_neighSolution[" << i << "][" << j << "] = " << ContTabu.m_neighSolution[i][j] << endl;
    }
  }

  return ( stream );
}
